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A  three-dimensional  model  of  polymer  electrolyte  fuel  cells  (PEFCs)  is  developed  to  investigate  multiphase 
flows,  species  transport,  and  electrochemical  processes  in  fuel  cells  and  their  interactions.  This  two-phase 
model  consists  of  conservation  principles  of  mass,  momentum,  species  concentration  and  charges,  and 
elucidates  the  key  physicochemical  mechanisms  in  the  constituent  components  of  PEFCs  that  govern  cell 
performance.  Efforts  are  made  to  formulate  two-phase  transport  in  the  anode  diffusion  media  and  its 
coupling  with  cathode  flooding  as  well  as  the  interaction  between  single-  and  two-phase  flows.  Numer¬ 
ical  simulations  are  carried  out  to  investigate  multiphase  flow,  electrochemical  activity,  and  transport 
phenomena  and  the  intrinsic  couplings  of  these  processes  inside  a  fuel  cell  at  low  humidity.  The  results 
indicate  that  multiphase  flows  may  exist  in  both  anode  and  cathode  diffusion  media  at  low-humidity 
operation,  and  two-phase  flow  emerges  near  the  outlet  for  co-flow  configuration  while  is  present  in  the 
middle  of  the  fuel  cell  for  counter-flow  one.  The  validated  numerical  tools  can  be  applied  to  investigate 
vital  issues  related  to  anode  performance  and  degradation  arising  from  flooding  for  PEFCs. 

Published  by  Elsevier  B.V. 


1.  Introduction 

Fuel  cells,  converting  the  chemical  energy  stored  in  fuels  directly 
and  efficiently  to  electricity  via  electrochemical  reactions,  have 
become  the  focus  of  new  energy  development  due  to  their  note¬ 
worthy  features  of  high  efficiency  and  low  emissions  [1-3].  Among 
all  types  of  fuel  cells,  the  polymer  electrolyte  fuel  cells  (PEFCs), 
also  called  polymer  electrolyte  membrane  (PEM)  fuel  cells,  have 
captured  the  public  attention  for  both  mobile  and  portable  appli¬ 
cations  [4,5].  In  addition  to  providing  high  power  density,  PEFCs 
work  at  low  temperature  (typically  <100  °C),  only  produce  water 
as  byproduct  and  can  be  compactly  assembled.  These  factors  make 
them  one  of  the  premiere  candidates  as  the  next  generation  power 
generator. 

Fig.  1  schematically  shows  a  single  PEFC.  A  typical  PEFC  consists 
of  the  following  components:  bipolar  plates  with  flow  channels 
grooved  in,  gas  diffusion  media  (DMs),  and  a  proton-conductive 
membrane  with  platinum  catalyst  coated  on  each  side.  The  bipo¬ 
lar  plates,  usually  graphite  or  metal  plates,  play  important  roles 
of  electronically  connecting  the  adjacent  cells  and  distributing  the 
reactant  gases  over  the  anode  and  cathode.  Gas  diffusion  media 
are  usually  carbon-based  substrates  such  as  carbon  paper  and 
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cloth  [6,7],  see  Fig.  2,  which  are  placed  on  the  membrane  elec¬ 
trode  assembly  (MEA)  and  perform  multiple  functions:  passages 
for  species  and  heat  transport  and  electronic  connection  between 
the  bipolar  plates  and  catalyst  layers.  Important  constituent  phases 
of  the  catalyst  layer  include  C/Pt  or  Pt  alloy  catalyst,  ionomer,  and 
void  space.  The  electrochemical  reactions  take  place  at  the  triple¬ 
phase  boundary  (TPB)  of  the  catalyst  layer.  The  proton-conducting 
membrane,  typically  made  of  Nation®,  plays  dual  roles  as  the  gas 
separator  and  electrolyte.  The  molecule  of  Nation®  is  characterized 
by  hydrophobic,  fluorinated  main  chains  with  hydrophilic  sulfonic 
acid  groups,  allowing  protons  to  weakly  attract  to  the  S03~  and 
travel  in  hydrated  region. 

Mathematical  modeling  of  PEFCs  has  been  a  rapidly  growing 
field  of  research.  Early  models  mostly  focus  on  electrochemi¬ 
cal  modeling  in  one  or  pseudo-two  dimensions  and  deal  with 
single-phase  phenomena  in  fuel  cells.  Later,  single-phase  multi¬ 
dimensional  models  were  developed,  e.g.  Garau  et  al.  [8],  Dutta 
et  al.  [9,10],  Um  et  al.  [11],  Mazumder  and  Cole  [12],  Meng  and 
Wang  [13],  and  Wang  and  Wang  [14].  Their  work  is  primarily 
based  on  the  computational  fluid  dynamics  (CFD)  approach,  where 
multi-dimensional  solutions  were  obtained  by  solving  transport 
equations  governing  conservation  of  mass,  momentum,  species, 
and  charge.  One  major  assumption  made  in  several  models  is  to 
treat  the  catalyst  layers  as  an  interface  without  thickness.  This 
assumption  simplifies  fuel  cell  modeling  and  facilitates  numerical 
implementation,  however  may  introduce  inaccuracy  due  to  neglect 
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Fig.  1.  Schematic  of  a  PEFC. 


of  electrolyte  phase  potential  variation  in  the  through-plane  direc¬ 
tion  of  the  catalyst  layer.  A  challenge  in  fuel  cell  modeling  is 
multiphase  transport  in  fuel  cells.  Multiphase  flow  originates  from 
water  production  by  the  oxygen  reduction  reaction.  Liquid  water 
affects  gaseous  reactant  supply  and  electrochemical  catalyst  activ¬ 
ity.  Two-phase  modeling  have  been  attempted  by  He  et  al.  [15], 
Wang  et  al.  [16],  Janssen  [17],  You  and  Liu  [18],  Berning  and  Djilali 
[19],  Pasaogullari  and  Wang  [20,21  ]  and  Weber  and  Newman  [22], 
which  mostly  dealt  with  isothermal  conditions.  Incorporation  of 
heat  transfer  has  been  attempted  by  Costamagna  [23],  Mazumder 
and  Cole  [24],  Birgersson  et  al.  [25],  Wang  and  Wang  [26],  Weber 
and  Newman  [27],  and  Meng  [28].  These  previous  models  mostly 
focus  on  steady-state  phenomena.  Dynamic  responses  of  fuel  cells 
have  been  investigated  by  Wang  and  Wang  [29],  Hu  and  Fan  [30], 
Shah  et  al.  [31  ],  Meng  [32],  and  Wang  [33]. 

Though  great  efforts  have  been  made  in  PEFC  modeling,  com¬ 
prehensive  descriptions  of  fuel  cells  that  elucidate  the  complex 
couplings  of  electrochemical  kinetics  and  multiphase  transport  are 
still  highly  in  need.  In  addition,  the  majority  of  previous  modeling 
work  only  address  two-phase  flow  at  high-humidity  operations, 
while  the  physics  of  transition  from  single  phase  to  multiple 
phases  and  formation  of  the  evaporation/condensation  fronts,  are 
extremely  valuable  and  of  paramount  importance  for  fuel  cell 
operation,  particularly  at  low  humidity.  Further,  most  of  previ¬ 
ous  studies  focus  on  the  cathode  flooding,  while  few  include  the 
two-phase  flow  in  the  anode.  Note  that  the  anode  electrode  may 


contribute  a  major  portion  of  the  ionic  resistance  [29],  particularly 
at  low  humidity. 

In  this  paper,  we  seek  to  develop  such  a  two-phase  fuel 
cell  model  to  elucidate  two-phase  transports  in  both  anode  and 
cathode  diffusion  media  and  their  interactions  as  well  as  the 
interactions  between  single-  and  two-phase  flows  and  between 
the  two-phase  flow  and  electrochemical  reaction.  3D  numerical 
simulations  are  carried  out  for  a  single  straight-channel  PEFC  at 
low-humidity  operation.  In  particular,  we  will  investigate  transi¬ 
tion  of  single-phase  flow  to  multiphase  one  and  formation  of  the 
evaporation/condensation  front  for  both  co-  and  counter-flow  con¬ 
figurations.  Simulation  results  are  presented  to  reveal  distributions 
of  two-phase  flow  in  various  parts  of  the  PEFC  and  validations  is 
made  for  both  low-  and  high-humidity  conditions. 

2.  Mathematical  model 


Fig.  1  schematically  shows  the  geometry  of  a  PEFC  and  the  con¬ 
stituent  components  to  be  modeled  in  this  work.  The  two-phase 
model  considers  the  electrochemical  and  transport  mechanisms 
in  all  of  the  key  components  including  the  catalyst  layers  and 
membrane.  The  following  assumptions  are  made :  ( 1 )  ideal  gas  mix¬ 
tures;  (2)  isotropic  and  homogeneous  membrane,  catalyst  layers 
and  gas  diffusion  media;  (3)  laminar  flow  due  to  small  pressure 
gradients  and  flow  velocity;  (4)  equilibrium  between  ionomer 
and  surrounding  fluid  in  the  catalyst  layers;  (5)  isothermal  con¬ 
dition.  The  model  consists  of  four  principles  of  conservation:  mass, 
momentum,  species  and  charge,  and  can  be  presented  in  concise 
form  as  follows: 


V  •  r  =  S,  where  r 

up 

uup  -  x 
YkUCk+Gk,di  ff 

(mfi)  C(g)V  .  n~ 

ywuCw  +  Gw?diff+  f  I  /^+Cw,perm+-^-^m^ 

-or^ffV0(m) 

-<TseffV0(s) 


and  5  = 


5c,< 

5CW 


S<£(m) 

-  %s)  _ 


(1) 


where  p,  u,  p,  Q,  Cw,  and  @^s\  denote  the  multiphase 

mixture  density,  superficial  fluid  velocity  vector,  pressure,  molar 
concentration  of  reactant  k,  water  molar  concentration,  electronic 


Fig.  2.  SEM  micrographs  of  (a)  carbon  paper  and  (b)  carbon  cloth  [7]. 
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Table  1 

Source  terms  of  the  conservation  equations  in  different  fuel  cell  regions 


phase  potential,  and  electrolyte  phase  potential.  Gdiff  represents  the 
species  diffusion  flux  in  gaseous,  liquid  and  solid  electrolyte  phases, 
while  Gw, perm  denotes  the  water  flux  due  to  hydraulic  permeation. 
The  electrochemical  and  transport  processes  are  coupled  together 
through  model  parameters  such  as  diffusion  coefficients  and  the 
source  terms,  Su,  SCk,  SCw,  S0( m),  and  S0(s ).  The  expression  of  these 
source  terms  is  summarized  in  Table  1  in  detail.  For  modeling  pur¬ 
pose,  the  pressure  gradient  term  in  the  momentum  equation  of  Eq. 
(1 )  is  included  in  the  source  term,  Su,  to  explicitly  show  the  expres¬ 
sion  of  the  Darcy’s  law  in  the  porous  diffusion  media  and  catalyst 
layers.  Discussion  of  these  source  terms  and  necessary  constitutive 
relations  as  well  as  the  key  electrochemical/transport  mechanisms 
are  elaborated  below. 

2.1.  Electrochemical  kinetics 

The  catalyst  layer  consists  of  four  components:  electrocatalyst, 
ionomer,  carbon,  and  void  space.  Platinum  and  the  alloys  of  plat¬ 
inum  and  ruthenium  are  the  typical  catalyst  materials  for  PEFCs. 
Currently,  Pt  loading  is  around  0.45  mg  cm-2  and  other  novel  tech¬ 
niques  have  been  reported  in  the  literature  to  reduce  the  value  to 
<0.2  mg  cm-2  [34,35].  The  following  electrochemical  reactions  take 
place  at  the  triple-phase  boundary  of  the  catalyst  layer: 

Flydrogen  oxidation  reaction  (FIOR)  in  the  anode: 

H2^  2H+  +  2e“  (2) 

Oxygen  reduction  reaction  (ORR)  in  the  cathode: 

02  +4e_  +4H+  2H20  (3) 


approximated  by  Tafel  kinetics: 

/  \  1/2 

In  the  anode  :  ja  =  af0e[  f  J  (^p^) 

In  the  cathode  :  jc  =  -a?0e[  exp  (~W v) 

where  the  surface  overpotential  is  defined  as 
r,  =  <J>(s)  -  <Z>(m)  -  Uo 


(6) 

(7) 

(8) 


The  equilibrium  potential,  U0  (V),  is  zero  in  the  anode,  while  in  the 
cathode  it  is  a  function  of  temperature  [36]: 


U0  =  1.23-  0.0009 (T  -  298)  (9) 

The  electrolyte  phase  potential,  <p(m\  is  the  driven  force  for 
proton  movement.  The  proton  conductivity,  crm,  is  determined 
by  local  membrane  water  content.  In  molecular  scale,  the  poly¬ 
mer  electrolyte  membrane,  typically  Nation®,  is  characterized  by 
hydrophobic,  fluorinated  main  chains  with  hydrophilic  sulfonic 
acid  side  chains.  The  protons  are  weakly  attracted  to  the  S03  -  group 
in  the  hydrated  region  and  are  able  to  travel  in  the  solid  electrolyte. 
The  relationship  between  the  ionic  conductivity  and  water  content 
in  a  Nation®  membrane  was  experimentally  measured  by  Springer 
etal.  [37]: 


crm  =  (0.  5139A-  0.326)  exp 


,26s  (3® 


(10) 


where  X  denotes  the  water  content  in  the  membrane,  defined  as  the 
ratio  of  the  number  of  water  molecules  to  the  number  of  charge 
(S03_H+)  sites.  We  adopt  the  following  formula  to  calculate  the 
value  of  A  [37]: 


f  0.043  +  17.81a  -  39.85a2  +  36.0a3  forO  <  a  <  1 
14  +  8s  forO  <  s  <  1  ’ 


Gw  ,  r  (T v  Psat(T) 

—  and  csat(T)  =  — 

where  log10Psat  =  -2.1794  +  0.02953(T  -  273.15)  -  9.1837 
x  10“5(T  -  273. 15)2  +  1.4454 
x  10“7(T-273.15)3  (11) 


The  electrochemical  kinetics  of  the  above  two  reactions  are 
described  by  the  well-known  Butler-Volmer  equation: 

j  =  ai0  | exp  (^Fr/j  -  exp  (-||Fp)  )  (4) 

where  aa  and  ac  are  anodic  and  cathodic  charge  transfer  coeffi¬ 
cients,  respectively.  The  exchange  current  density,  z0,  is  determined 
by  the  catalyst  electrochemical  kinetics,  while  the  surface-to- 
volume  ratio,  a,  describes  the  roughness  of  porous  electrodes.  The 
presence  of  liquid  water  in  the  catalyst  layer  may  reduce  the  reac¬ 
tion  area  and  the  following  empirical  formula  is  generally  adopted 
to  account  for  the  liquid  water  coverage: 

a  =  (l-s)Tca0  (5) 

in  which  the  liquid  water  saturation,  s,  is  defined  as  volume  fraction 
of  liquid  water  in  the  void  space. 

In  PEFCs,  FIOR  is  fast,  thus  yielding  a  low  anode  overpotential. 
Therefore  Eq.  (4)  in  the  anode  can  be  adequately  simplified  to  a 
linear  kinetic  equation.  For  ORR,  sluggish  kinetics  results  in  a  high 
cathode  overpotential.  Thus,  the  Butler-Volmer  equation  can  be 


where  a  denotes  the  water  vapor  activity.  In  catalyst  layers  or 
Gore™  membranes,  the  value  of  the  ionic  conduction  coefficient  is 
modified  according  to  the  Bruggeman  relation: 

°'mf  =  Sm°'m  (12) 

where  rm  is  the  Bruggeman  factor  and  £m  denotes  the  volume  frac¬ 
tion  of  ionomer.  In  addition,  the  work  of  Ma  et  al.  [38]  indicates 
ionic  conductivity  anisotropy  in  the  membrane.  In  that  case, 
in  Eqs.  (10)  and  (12)  can  be  written  in  a  vector  form  denoting  the 
conduction  coefficient  in  three  dimensions. 

Once  the  electrolyte  phase  potential  is  calculated  through  Eq. 
(1 ),  the  protonic  current  flux  in  the  membrane  can  be  expressed  as 

i<m)  =  -a8ffV$(m>  (13) 

Note  that  i-m)  is  a  vector  and  represents  the  current  flows  in  three 
dimensions.  To  calculate  the  average  current  density,  one  can  inte¬ 
grate  z(m)  over  a  cross-section  of  the  membrane: 

1=2-  JJ  ?<m>-dS  (14) 

cross-section 
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2.2.  Multiphase  flow 

One  of  the  most  complex  phenomena  in  PEFCs  is  multiphase 
transport,  originating  from  water  production  by  the  oxygen  reduc¬ 
tion  reaction.  Liquid  droplets  may  block  pore  paths  and  hamper  the 
reactant  transport  to  the  reaction  site,  leading  to  a  substantial  volt¬ 
age  loss.  Following  the  M2  theory  [39],  we  define  the  two-phase 
mixture  density  as 

p  =  sp(1)  +  (1  -  s)p(g)  (15) 


Gdiff  =  -D(s>’effVC(s),  where  D^’eff  = -D^  =  and 

Note  that  the  term  on  the  right  side  of  the  second  equation  is  also 
referred  to  as  the  Bruggeman  relation.  In  the  multiphase  region,  the 
liquid  attaches  on  the  pores’  wall,  following  the  same  morphology 
of  the  DM  solid  matrix.  Thus,  the  effective  gas  diffusion  coefficient 
is  modified  by 


where  s  is  the  liquid  water  saturation.  In  addition,  p(g)  is  determined 
by  the  constituent  species  and  their  contents  in  the  gas,  and  its 
values  on  the  anode  and  cathode  sides  are  quite  different.  Note  that 
the  mixture  density  is  the  gaseous  density  when  no  liquid  water 
is  present.  The  liquid  saturation,  s,  is  obtained  from  the  following 
relation  with  the  mixture  water  concentration,  Cw,  after  the  latter 
is  solved  from  the  water  equation  in  Eq.  (1): 


s  = 


0 

Cw  —  CSat 
p(')/Mw  -  Csat 


Cw  <  Qat 


Cw  >  CSat 


(16) 


D(s)>eff  =  [g(i  -  s)]TdD(s) 


(19) 


Liquid  flows  also  affect  species  transport  and  the  impact  is 
included  via  the  convection  corrector  factor,  y. 


M(g) 

Yk  ~  p(s)(l  -s) 


X«  xfe)  \ 

Mw  +  ~^  ) 


(20) 


where  the  relative  mobilities  of  individual  phases,  are 


The  interaction  of  the  two  phase  flows  in  the  diffusion  media  is 
described  through  the  relative  permeabilities,  and  k[g\  defined 
as  the  ratio  of  the  intrinsic  permeability  of  liquid  and  gas  phases, 
respectively,  to  the  total  intrinsic  permeability  of  a  porous  medium. 
Physically,  these  parameters  describe  the  extent  to  which  one  fluid 
is  hindered  by  others  in  pore  spaces,  and  hence  can  be  formulated 
as  a  function  of  liquid  saturation.  Most  of  previous  work  adopted 
cubic  relations.  Here,  we  take  the  following  formula  for  relative 
permeabilities  [40]: 

k^  =  s4  and  k'g)  =  (l-s)4  (17) 


k®/v® 


(k™/vffl)  +  (k[.s)/v(s)) 


and  A/8*  = 


kW/v(s) 


(k™/v’(l))  +  (k‘s)/v(g)) 


(21) 


Liquid  water  transport  in  the  diffusion  media  is  driven  by  the  liquid 
pressure  which  is  calculated  by 


P(l)  =  Pc 


/  P  \  V  2 

+  P(g)  =  rcos(Pc)f  J{s)  +  P(g) 


(22) 


In  pure  gas  region,  s  is  equal  to  zero  therefore  the  above  multiphase 
flow  equation  is  back  to  the  single-phase  expression. 

In  addition,  the  flows  may  also  be  affected  by  the  electrochemi¬ 
cal  activities  in  the  catalyst  layers  as  well  as  the  transport  processes 
between  the  solid  electrolyte  and  fluid.  A  mass  source  term  is 
generally  added  in  Eq.  (1)  to  describe  this  effect  as  derived  by 
Wang  and  Wang  [14],  who  also  concluded  that  the  mass  source 
term  has  negligible  impacts  on  reactant  transport  and  cell  per¬ 
formance  prediction.  Similar  assumption  of  neglecting  this  source 
term  has  been  adopted  by  both  single-  and  two-phase  models  of 
fuel  cells  [7,11,12,20,21,29]  with  some  of  them  validated  experimen¬ 
tally  [7,11].  A  recent  paper  has  presented  a  theoretical  discussion 
on  this  assumption  [41  ]  and  raised  some  concerns  regarding  its 
validity. 

2.3.  Reactant  and  water  transport 

In  gaseous  phase,  the  species  of  hydrogen,  oxygen  and  water 
follow  similar  transport  mechanisms.  For  the  traditional  design  of 
PEFC  flow  fields,  the  convective  mass  transfer  dominates  the  species 
transport  in  gas  channels,  while  diffusion  is  the  major  mechanism 
in  porous  media.  For  the  interdigitated  flow  field  [42]  or  the  one  in 
Ref.  [43],  convection  can  be  the  dominant  mechanism  in  diffusion 
media. 

2.3.1.  Transport  in  gaseous/liquid  phase 

The  multi-component  diffusion  in  the  gaseous  phase  is 
described  by  the  Stefan-Maxwell  equation.  To  simplify  the  expres¬ 
sion,  the  Fields  law  is  generally  used.  The  diffusion  flux,  Gdiff,  and 
the  associated  diffusion  coefficient  in  various  components  of  fuel 
cells  can  be  unified  as 


where  Pc  is  the  capillary  pressure  and  r  is  the  surface  tension.  The 
expression  of  Pc  in  the  above  is  called  the  Leverett  function,  and  J{s) 
for  hydrophobic  diffusion  media  is  given  by 

J(s)  =  1.417s  -  2.120s2  +  1.263s3  (23) 

Note  that  the  Leverett  /-function  only  considers  the  impact  of  two 
characteristics  of  a  porous  medium,  i.e.  porosity  and  permeabil¬ 
ity,  while  ignoring  the  effect  of  detailed  pore  morphology  [44].  In 
addition,  diffusion  media  materials  are  typically  made  hydrophobic 
through  adding  PTFE  to  facilitate  water  removal.  The  PTFE  load¬ 
ing,  commonly  ranging  from  5  to  30%,  has  significant  influence 
on  the  contact  angle,  0C.  High  PTFE  loadings  have  been  studied 
by  Lin  and  Nguyen  [45].  Once  the  capillary  pressure  is  calculated, 
the  flux,  p\  in  the  water  equation  of  Eq.  (1)  can  be  obtained 
through 

p  =  ^fK[VPc  +  (p(l)  _  p(g))g]  (24) 

Note  that  the  capillary  flux  term  ((mfj^/Mw)  -  (c[s)/p^g))}/^  in  the 
water  equation  of  Eq.  (1)  also  accounts  for  the  water  transport  by 
the  gas  flow  induced  by  capillary  liquid  flow  [39]. 

In  addition,  note  that  the  liquid  water  transport  in  the  chan¬ 
nel  may  affect  the  two-phase  transport  in  the  diffusion  media 
through  the  interaction  at  the  interface  between  the  diffusion 
media  and  channel.  However,  at  current,  few  two-phase  models 
have  been  developed  to  accurately  describe  the  channel  flood¬ 
ing  due  to  its  complexity.  Therefore,  most  of  previous  models 
[7,20,21,24]  including  the  one  in  this  paper  neglect  its  impacts 
and  assume  zero  liquid  water  saturation  on  diffusion  media  sur¬ 
face,  with  focus  on  the  two-phase  physics  within  the  diffusion 
media. 
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2.3.2.  Water  transport  in  the  ionomer 

Water  transport  also  takes  place  in  the  solid  electrolyte.  Water 
in  the  membrane  is  essential  for  the  ionic  conductance,  see  Eq.  (10). 
Two  processes  may  affect  the  membrane  hydration,  one  is  the  water 
electro-osmotic  drag;  the  other  is  water  back-diffusion.  The  coef¬ 
ficient  of  the  electro-osmotic  drag,  nd,  depends  on  the  local  water 
content  [46]: 


for  A  <  14 
(A  -  14)  +  1 .0  otherwise 


(25) 


The  water  diffusion  coefficient  in  the  membrane  is  also  a  function 
of  water  content  [47]  and  the  diffusive  flux  is  given  by 

^Wjdiff  —  v  Lw 


where 


D(m)  _  f  3.1 
W  f  4.1 


x  10  3A(e°-28A  -  i)e-2436/T  for  0  <  A  <  3 


17  x  10  4A(1  +  161  e-A)e_2436/T  otherwise 


ficient  in  the  portion  of  the  ionomer  is  modified  by 

D^}’eff  =  (28) 

Assuming  local  microscopic  thermodynamic  equilibrium 
between  ionomer  and  surrounding  fluid,  one  can  combine  the 
diffusive  transport  in  the  void  space  and  ionomer  by  defining  an 
effective  diffusion  coefficient: 

£)eff  _  g^d  D (g)  -L  sr m  ^  ^  RT  dx  (m) 

uw  -£  uw  +sm  EW  p ^  dQUw  (2 yj 

where  is  the  density  of  a  dry  membrane.  Note  that  by  defining 
thermodynamic  equilibrium  we  combine  water  transport  equa¬ 
tions  in  different  phases  and  there  is  no  need  to  explicitly  express 
the  water  exchange  rate  at  the  interface  in  the  final  form  of  the 
equation. 

In  addition  to  diffusion,  liquid  water  pressure  difference 
between  the  two  membrane  surfaces  may  lead  to  hydraulic  perme¬ 
ation  through  the  membrane.  The  permeation  flux  is  determined 
by  the  permeability  of  the  membrane,  I<m,  and  liquid  pressure  gra¬ 
dient: 

*  Krr 


MwvC) 


vp(') 


(30) 


Here,  we  assume  liquid  water  pressure  in  the  membrane  is  lin¬ 
ear  in  the  thickness  direction,  determined  by  the  pressures  at  the 
membrane  surface. 

3.  Boundary  conditions 

Eq.  (1)  form  a  complete  set  of  governing  equations  with  eight 
unknowns:  u  (three  components),  P,  Q,  Cw,  <P(m),  and  Their 
corresponding  boundary  conditions  are  described  as  follows. 

3.1.  Flow  inlet  boundaries 

The  inlet  velocity  uin  in  a  gas  channel  is  expressed  by  the  respec¬ 
tive  stoichiometric  flow  ratio,  i.e.  §a  or  fc.  defined  at  the  average 
current  density,  /,  as 


(  la  \ 

/Am 

.  ,  “  F“ 

inlet 

2r 

V  4CoA  / 

(31) 


inlet 


where  Aa,  Ac,  and  Am  are  the  flow  cross-sectional  areas  of  the  anode 
and  cathode  gas  channels  and  the  membrane  area,  respectively.  The 


inlet  molar  concentrations  are  determined  by  the  inlet  pressure  and 
humidity  according  to  the  ideal  gas  law. 

3.2.  Outlet  boundaries 

Fully  developed  or  no-flux  conditions  are  applied. 

3.2.1.  Walls 

No-slip  and  impermeable  velocity  condition  and  no-flux  condi¬ 
tion  are  applied: 


; ' 
Q 


(26) 

3  n 

Cw 

0(m) 

\  0IS>  J 

(27) 

(  p  \ 

3 

Q 

coef- 

3  n 

Cw  1 

^0(m)  J 

=  0,  u  I  wail  =  0,  and 


outlet 


=  o 


(32) 


wall 


The  boundary  conditions  for  the  electronic  phase  potential,  &(s\  at 
the  bipolar  plate  outer  surfaces  can  be  expressed  as 

30(s) 


^^lanode  —  0, 


3  n 


"  I  cathode  — 


/An 

reffA 


and 


c.wall 


3  n 


"  I  otherwise  —  0 


(33) 


where  Acwall  is  the  area  of  the  cathode  outer  surface. 

4.  Numerical  procedures 

The  governing  equation,  Eq.  (1),  along  with  their  appropriate 
boundary  conditions  is  discretized  by  the  finite-volume  method, 
with  SIMPLE  (semi-implicit  pressure  linked  equation)  algorithm 
[48].  For  finite-volume  discretization,  it  is  convenient  to  unify  all 
governing  equations,  including  the  transient  terms,  in  the  following 
form: 


vr((9)  =  S© 


(34) 


where  0  stands  for  any  dependent  variable  in  Eq.  (1).  Integrating 
the  above  equation  throughout  an  arbitrary  volume  V  bounded  by 
a  closed  surface  S,  yields: 


r(0)dS=  /  S@dv 

Jv 


(35) 


where  S  is  the  surface  vector.  Taking  V  and  S  to  be  the  volume  Vp 
and  discrete  faces,  Sj  t  of  a  computational  cell,  respectively,  one  can 
reach: 


£  f  H&)  ■  dS  =  [  S0  dv 
.  Js,  Jv, 


(36) 


J  Jsi 


The  final  form  of  the  discrete  finite-volume  equation  can  be 
expressed  as 


Bp0g  =  £Bm0£+B(0°) 


(37) 


The  above  equation  is  then  solved  by  the  algebraic  multi-grid  (AMG) 
method. 

The  mesh  of  a  single-channel  PEFC  employed  here  for  a  numer¬ 
ical  study  is  shown  in  Fig.  3.  About  120,000  (60  x  100  x  12) 
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Fig.  3.  Computational  domain  and  mesh  of  a  single-channel  PEFC. 


Table  2 

Geometrical,  physical,  and  operating  parameters 


Quantity 

Value 

Gas  channel  depth/width 

0.5/1.0mm 

Shoulder  width 

1.0  mm 

GDL/catalyst  layer/membrane  thickness,  8 

0.2/0.01/0.03  mm 

Anode/cathode  pressures,  P 

2.0/2.0atm 

Stoichiometry,  £a/£c  @  0.8  A  cm-2 

1.5/2.0 

Fuel  cell  temperature,  T 

353.15  K 

Porosity  of  GDLs  [33]/catalyst  layers,  e  [40] 

0.6/0.4 

Volume  fraction  of  ionomer  in  catalyst  layers,  em 

0.2 

[7] 

Tortuosity,  rd/rm  [11] 

1.5/1. 5 

Electronic  conductivity  of  GDLs/bipolar  plates,  crff 

500/2000  Wm-1  K-1 

[29] 

Viscosity  of  liquid  water,  /zO)  [20,26] 

3.5  x  10-4  kgm-1  s-1 

Surface  tension,  liquid-water-air  (80  °C),  r  [20,26] 

0.0625  Nm-1 

Contact  angle,  0C  [7] 

110° 

Permeability  of  GDL,  I<DMII<m  [26] 

10_12/5  x  10-20  m2 

Exchange  current  density  x  reaction  surface  area, 

1.0  x  109/3.5  x  104  Am-3 

ao*'o,a/aoio,c  [7] 

Species  diffusivity  in  anode  gas  @standard 

U/U  x  lO"4  m2  s-1 

condition,  D0  //  /w  [29] 

Species  diffusivity  in  cathode  gas  @standard 

3.24/3.89  x  10-5  m2  s"1 

condition,  D0  o2/w  [29] 

computational  cells  are  used  to  capture  the  complex  electrochemi¬ 
cal  and  physical  phenomena  in  the  PEFC.  Geometrical  and  operating 
parameters  of  the  PEFC  as  well  as  physical  properties  are  listed 
in  Table  2.  An  average  current  density  is  specified  as  an  input 
parameter,  allowing  the  local  current  density  and  electronic  phase 
potential  to  vary  spatially  according  to  local  conditions.  In  all  the 
simulations  to  be  presented  in  the  next  section,  values  of  equation 
residuals  are  smaller  than  10-6. 

5.  Results  and  discussion 

A  single-channel  PEFC  with  the  Gore™  30  p,m  membrane  and 
carbon  cloth  diffusion  media  is  chosen  for  a  case  study.  Both  co- 
and  counter-flow  configurations  of  the  anode  and  cathode  streams 
are  considered.  Air  and  pure  hydrogen  are  fed  in  the  PEFC.  All 
results  are  intended  to  reveal  and  explore  the  two-phase  phenom¬ 
ena  in  a  PEFC  operating  at  0.8  A  cm-2  and  inlet  humidification  of 
RHa/c  =  66/66%.  This  operation  condition  indicates  typical  patterns 
of  two-phase  phenomena  in  both  anode  and  cathode,  which  are 
sufficient  to  demonstrate  the  fundamentals  of  two-phase  trans¬ 
port  within  the  diffusion  media  as  well  as  the  interaction  between 


anode  and  cathode  two-phase  flows  in  both  co-  and  counter-flow 
configurations. 

Fig.  4  shows  contours  of  the  liquid  water  saturation  in  the  diffu¬ 
sion  media  under  both  channel  and  land.  It  can  be  seen  that  single- 
and  multiphase  regions  coexist  in  PEFCs  at  this  low-humidity  oper¬ 
ation.  The  single-phase  region  is  near  the  inlet  where  the  dry  gases 
are  fed  in.  Due  to  water  production  by  fuel  cells,  liquid  water 
emerges  downstream  and  the  flow  in  the  diffusion  media  shifts 
to  gas-water  multiphase  flow.  Liquid  water  emerges  first  in  the 
cathode  side  and  higher  saturation  levels  appear  under  the  cath¬ 
ode  land  with  the  value  as  high  as  ~25%.  In  addition,  despite 
the  water  electro-osmotic  drag,  there  exists  slight  flooding  in  the 
anode,  which  may  affect  the  anode  electrochemical  activity  and 
reactant  transport.  Anode  flooding  was  also  shown  in  the  exper¬ 
iment  of  Ref.  [49].  Further,  with  the  increase  of  the  flow  rate  in 
the  channel  or  stoichiometric  ratio,  the  channel  stream  can  remove 
more  water  and  hence  reduce  the  two-phase  region  area  till  liquid 
disappears.  In  that  case,  the  present  two-phase  model  changes  back 
to  single-phase  one  and  such  single-phase  phenomena  has  been 
addressed  adequately  in  previous  literature  [12,14,43].  In  addition, 
lowering  the  external  humidification  will  also  diminish  the  two- 
phase  region. 

Fig.  5  shows  the  reactant  concentration  contours  in  the  mid¬ 
plane  of  gas  channels  and  diffusion  media.  It  can  be  seen  that 
oxygen  and  hydrogen  concentrations  decrease  down  the  channel 


Fig.  4.  Liquid  water  saturation  distributions  in  the  diffusion  media:  (a)  under  the 
channel  and  (b)  under  the  land  (operating  conditions:  Pa/C  =  2  atm,  Tceu  =  80°C, 
RHa/c  =  66/66%,  and  Stoich.a/c  =  1.5/2). 
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Fig.  5.  Reactant  profiles  of  (a)  hydrogen  and  oxygen  contours  at  the  middle  section  of  the  PEFC  and  (b)  hydrogen  and  oxygen  concentrations  under  the  land  and  channels  at 
the  A-A  location  (operating  conditions:  Pa/C  =  2  atm,  Tcell  =  80  °C,  RHa/c  =  66/66%,  and  Stoich.a/c  =  1.5/2). 


tLx  fjo- 


Fig.  6.  Relative  humidity  (RH)  distributions:  (a)  at  the  middle  section  of  the  fuel  cell  under  the  channel;  (b)  at  the  cross-section  of  the  plane  A-A  and  (c)  at  the  cross-section 
of  the  plane  B-B.  The  gray  areas  denote  the  multiphase  region.  Here  RH  is  defined  by  the  same  formula  as  the  water  activity,  a  (operating  conditions:  Pa/C  =  2  atm,  Tce\\  =  80  °C, 
RHa/c  =  66/66%,  and  Stoich.a/c  =  1.5/2). 
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Fig.  7.  The  velocity  and  saturation  distributions  in  the  cathode  GDL:  (a)  liquid  water  velocity  and  saturation  at  the  cross-section  of  50%  fraction  distance  from  the  inlet;  (b) 
gas  velocity  at  the  cross-section  of  50%  fraction  distance  from  the  inlet;  (c)  liquid  water  velocity  and  saturation  at  the  cross-section  of  90%  fraction  distance  from  the  inlet; 
(d)  gas  velocity  at  the  cross-section  of  90%  fraction  distance  from  the  inlet  (operating  conditions:  Pa/C  =  2  atm,  Tcell  =  80  °C,  RHa/c  =  66/66%,  and  Stoich.a/c  =  1.5/2). 


due  to  reaction  consumption.  In  contract  to  the  substantial  decrease 
in  the  along-channel  direction,  the  concentration  only  experiences 
a  small  decline  across  the  GDLs.  In  addition,  the  transport  resistance 
under  the  land  is  relatively  large,  leading  to  a  considerable  drop 


in  the  reactant  concentration.  As  the  hydrogen  diffusivity  is  larger 
than  the  one  of  oxygen  in  the  air,  the  hydrogen  profile  is  more  even 
in  the  through-plane  direction  and  its  concentration  undergoes  a 
smaller  drop  under  the  land. 


Fig.  8.  (a)  Schematic  of  the  PEFC  with  counter-flow  configuration,  and  the  contours  of  liquid  water  saturation:  (b)  under  the  channel  and  (c)  under  the  land  (operating 
conditions:  Pa/C  =  2  atm,  Tcell  =  80  °C,  RHa/c  =  66/66%,  and  Stoich.a/c  =  1.5/2). 
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Fig.  9.  The  distributions  of  (a)  the  liquid  water  saturation,  (b)  liquid  water  velocity,  and  (c)  gaseous  velocity  at  middle  section  of  the  cathode  GDL  for  the  counter-flow 
(operating  conditions:  Pa/C  =  2  atm,  Tcen  =  80 °C,  RHa/c  =  66/66%,  and  Stoich.a/c  =  1.5/2). 


Fig.  6  shows  the  relative  humidity  (RH)  profiles  at  different  loca¬ 
tions  of  the  fuel  cell.  In  Fig.  6(a),  the  RH  in  the  anode  side  decreases 
initially  due  to  the  water  loss  caused  by  the  electro-osmotic  drag. 
It  starts  increasing  as  back-diffusion  dominates  water  transport 
across  the  membrane.  Here,  RH  is  defined  by  the  same  formula 
as  the  water  activity,  a,  in  Eq.  (11),  which  considers  the  water 
molecules  in  both  liquid  and  gaseous  phases,  therefore  the  RH  can 
be  over  1  as  shown  in  the  figure.  Fig.  6(b)  shows  the  RH  distribu¬ 
tion  at  the  cross-section  of  plane  A-A,  where  liquid  water  emerges. 
At  this  location,  the  anode  is  subjected  to  flooding  under  the  land 
while  there  is  no  liquid  water  under  the  channel.  At  the  location 
of  B-B,  i.e.  upstream  region,  the  RH  is  lower  than  1  on  both  side  as 
shown  in  Fig.  6(c)  and  therefore  no  liquid  water  appears.  In  contrast 
to  Fig.  6(b),  Fig.  6(c)  indicates  that  the  lowest  water  concentration 


in  the  anode  is  in  the  diffusion  media  under  the  channel,  therefore 
both  anode  channel  and  the  region  under  the  land  supply  water  to 
the  anode  catalyst  layer  under  the  channel  to  compensate  for  the 
water  loss  caused  by  the  electro-osmotic  drag.  This  is  partly  due  to 
the  higher  transport  resistance  under  the  land,  which  forces  more 
produced  water  in  the  cathode  to  diffuse  back  to  the  anode  side 
under  the  land. 

Fig.  7  shows  the  velocity  and  saturation  distributions  in  the 
cathode  GDL  at  two  cross-sections.  At  the  location  of  50%  fraction 
distance,  there  exists  slight  flood  under  the  channel,  therefore  the 
local  liquid  velocity  is  quite  small  as  shown  in  Fig.  7(a).  In  addi¬ 
tion,  the  liquid  flows  are  uniformly  from  the  catalyst  layer  to  the 
gas  channel  and  from  under-land  region  to  the  one  under  the  chan¬ 
nel,  which  shows  water  removal  from  the  catalyst  site.  The  mass 


Fig.  10.  The  contours  of  protonic  current  densities  at  the  middle  section  of  the  membrane  for  (a)  co-flow  and  (b)  counter-flow  cases  (operating  conditions:  Pa/C  =  2  atm, 
Tcen  =  80  °C,  RHa/c  =  66/66%,  and  Stoich.a/c  =  1.5/2). 
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Fig.  11.  Polarization  curves  under  both  low-  and  high-humidity  conditions.  The 
experimental  data  are  obtained  from  Ref.  [7]. 

flow  of  the  liquid  water  is  accompanied  by  a  reverse  gas  flow  in  the 
porous  media,  indicative  of  the  complexity  of  the  multiphase  flow 
phenomena  in  fuel  cells. 

Previous  results  are  for  the  operation  of  a  PEFC  with  co-flow  con¬ 
figuration.  Alternative  flow  field  currently  under  intensive  study  is 
counter-flow,  i.e.  the  fuel  and  air  flows  are  set  in  opposite  directions. 
Fig.  8(a)  shows  schematic  of  a  PEFC  with  such  flow  field.  The  idea 
of  the  counter-flow  design  is  to  humidify  the  inlet  reactant  flows 
through  an  internal  water  recirculation  and  therefore  reduce  the 
requirement  for  external  humidification.  Fig.  8(b)  and  (c)  present 
contours  of  liquid  water  saturation  in  the  middle  section  of  the 
under-channel  and  under-land  regions,  respectively.  In  contrast  to 
the  ones  of  co-flow  case  in  Fig.  4,  liquid  water  appears  only  in  the 
middle  of  the  fuel  cell.  Similarly,  the  cathode  is  subjected  to  a  more 
severe  flooding  and  liquid  water  first  emerges  in  the  cathode.  How¬ 
ever,  as  the  cathode  outlet  is  set  close  to  the  dry  anode  inlet  region, 
water  is  transported  to  the  anode  side  via  the  membrane,  leading 
to  the  transition  of  cathode  two-phase  flow  back  to  single-phase 
one.  As  the  anode  wet  flow  downstream  humidifies  the  dry  cath¬ 
ode  inlet  reactant,  liquid  water  emerges  much  earlier  in  the  cathode 
comparing  with  the  co-flow  case. 

Fig.  9  shows  distributions  of  liquid  water  saturation,  together 
with  velocities  of  liquid  and  gas  flows,  at  the  middle  section  of  the 
cathode  diffusion  media  for  the  counter-flow.  It  can  be  seen  that 
the  under-land  GDL  is  subjected  to  more  severe  flooding  and  liquid 
flow  is  uniformly  from  the  under-land  region  to  the  under-channel 
one  where  it  is  taken  away  by  the  channel  stream.  In  addition,  the 
channel  has  a  pressure  gradient  along  the  channel  due  to  the  chan¬ 
nel  viscous  flow,  which  will  affect  the  gaseous  flow  in  the  diffusion 
media  according  to  the  Darcy’s  law  and  induce  an  along-channel 
component  of  the  gas  velocity,  as  shown  in  Fig.  9(c). 

Fig.  10  shows  contours  of  the  current  densities  at  the  mid¬ 
dle  section  of  the  membrane  for  both  co-  and  counter-flow  cases. 
Both  cases  show  typical  distributions  of  the  current  density  under 
low-humidity  operation,  i.e.  higher  current  density  appears  in 
the  middle  of  the  fuel  cell.  Near  the  inlets,  dry  reactants  dehy¬ 
drate  the  membrane,  leading  to  a  substantial  ohmic  resistance. 
As  to  the  co-flow,  liquid  water  near  the  outlet  increases  the  mass 
transport  polarization,  leading  to  a  low  local  performance  down¬ 
stream.  Fig.  11  shows  the  validation  result  under  both  low-  and 
high-humidity  conditions.  The  output  voltage  is  a  key  parameter 
inductive  of  cell  performance.  The  figure  shows  a  good  agreement 
between  the  model  predictions  and  experimental  data. 

6.  Conclusions 

In  this  paper,  a  fuel  cell  model  was  developed  and  3D  numer¬ 
ical  simulations  were  carried  out  to  elucidate  the  fundamentals 


of  complex  two-phase  transport  in  both  anode  and  cathode  dif¬ 
fusion  media  and  their  interactions  as  well  as  the  interactions 
between  single-  and  two-phase  flows  and  between  two-phase 
flows  and  electrochemical  reactions.  Simulation  results  revealed 
detailed  two-phase  phenomena  under  low  humidity  operation  for 
both  co-  and  counter-flow  configurations.  We  found  that  two-phase 
flow  can  occur  in  both  anode  and  cathode  and  that  at  low-humidity 
operation  single-  and  multiphase  flows  coexist  in  a  fuel  cell.  In  the 
co-flow  configuration,  liquid  emerges  downstream  due  to  water 
production  by  fuel  cells,  while  flooding  is  more  severe  in  the  mid¬ 
dle  of  the  fuel  cell  for  the  counter-flow  configuration  due  to  the 
internal  humidification.  In  both  flow  configurations,  flooding  area 
under  the  land  is  larger  than  the  one  under  the  channel.  The  voltage 
predictions  of  the  model  agreed  well  with  experimental  data.  The 
validated  model  can  be  employed  for  detailed  fundamental  and 
parametric  studies,  such  as  key  parameters  governing  the  anode 
flooding,  impacts  of  the  anode  flooding,  factors  affecting  the  single- 
and  two-phase  regions,  and  advanced  water  management,  for  next 
generation  high  performance  PEFCs. 
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Appendix  A.  Nomenclature 


a  water  activity;  effective  catalyst  area  per  unit  volume 

(m2  m-3) 

a0  catalyst  surface  area  per  unit  volume  (m2  m-3 ) 

A  electrode  area  (m2) 

C  molar  concentration  of  species  k  (mol  m-3 ) 

D  species  diffusivity  (m2  s-1 ) 

EW  equivalent  weight  of  dry  membrane  (kg  mol-1 ) 

F  Faraday’s  constant  (96,487  C  per  equivalent) 

G  species  diffusion/permeation  flux  (mol  m-2 ) 

ie  superficial  current  density  (A  cm-2 ) 

I  current  density  (A  cm-2 ) 

j  transfer  current  density  (A  cm-3 ) 

j(1)  mass  flux  of  liquid  phase  (kg  m-2  s-1 ) 
kr  relative  permeability 

I<  permeability  (m2) 

L  length  (m) 

mass  fraction  of  species  k  in  liquid  phase 
M  molecular  weight  (kg  mol-1 ) 

n  the  direction  normal  to  the  surface 

nd  electro-osmotic  coefficient,  H20/H+ 

P  pressure  (Pa) 

R  universal  gas  constant  (8.134J  mol-1  K-1 ) 

s  liquid  saturation 

S  source  term 

t  time(s) 

T  temperature  (K) 

u  velocity  vector  (ms-1) 

U0  equilibrium  potential  (V) 

Greeks 

a  transfer  coefficient;  net  water  flux  per  proton  flux 

yc  correction  factor  for  species  convection 
8  thickness  (m) 

s  porosity 

r)  surface  overpotential  (V) 
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Oc  contact  angle  (°) 

0  phase  potential  (V) 

X  membrane  water  content 

X^  mobility  of  phase  k 

v  kinematic  viscosity  (m2  s-1 ) 

§  stoichiometric  flow  ratio 

p  density  (kg  m-3) 

a  electronic  or  ionic  conductivity  (S  m-1 ) 

r  surface  tension  (Nm-1) 

r  shear  stress  (Nm-2) 

Superscripts  and  subscripts 
a  anode 

c  cathode;  capillary 

CL  catalyst  layer 

D  diffusion 

DMs  diffusion  media 

eff  effective  value 

g  gas  phase 

GDL  gas  diffusion  layer 

in  inlet 

k  species;  liquid  or  gas  phase 

1  liquid 

m  membrane  phase 

o  gas  channel  inlet  value;  reference  value 

ref  reference  value 

s  solid 

sat  saturate  value 

w  water 
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